function W = integrate_W(dWdd,delta_grid,delta_grid_stp)

W = zeros(length(delta_grid),1);
% W(1) = 0;

h = delta_grid_stp;
for i = 2:length(delta_grid)
    W(i) = W(i-1) + ((dWdd(i) + dWdd(i-1))/2)*h;
end

end